1 Loading required packages

knitr::opts_chunk$set(warning = FALSE, message = FALSE)
library(tidyverse)
library(sf)
library(tmap)
library(PerformanceAnalytics) #scatter plot matrix
library(PerformanceAnalytics) #scatter plot matrix
 
# Change the presentation of decimal numbers to 4 and avoid scientific notation
options(prompt="R> ", digits=5, scipen=999)

2 Import the shapefile

sh<- st_read("../../Data/city.shp")
Reading layer `City' from data source `C:\Users\felip\OneDrive\Desktop\Meidai\Seminar Research\GitHub\book-2020-spatial-analysis-methods-and-practice\Data\City.shp' using driver `ESRI Shapefile'
Simple feature collection with 90 features and 15 fields
geometry type:  POLYGON
dimension:      XY
bbox:           xmin: 472360 ymin: 4199900 xmax: 481460 ymax: 4209300
projected CRS:  GGRS87 / Greek Grid
object.size(sh)
140128 bytes
s.sp <- as(sh, "Spatial")
class(s.sp)
[1] "SpatialPolygonsDataFrame"
attr(,"package")
[1] "sp"
class(sh)
[1] "sf"         "data.frame"

2.1 Inspect the table within the shapefile

st_drop_geometry(sh)

2.2 A quick map

qtm(sh)

3 Creating the breaks for the custom breaks map

breaks<- c(-Inf, 15000,20000,25000, Inf)

3.1 creating the custom breaks map

tm_shape(sh) +
  tm_fill("Income", title= "Income",breaks=breaks,palette="Greens")+
  tm_borders()

4 Creating a histogram

sh %>% 
  ggplot(aes(x=Income))+
  geom_histogram(bins=10, color="black", fill="white")

5 creating a box plot

geom_jitter(size = 3, alpha = 0.3, width = 0.1) is used so that the points are not all crammed in the same line

set.seed(2019)
sh %>% 
  ggplot(aes(x=Income, y=Municipali))+
  geom_boxplot()+
  geom_jitter(size = 3, alpha = 0.3, width = 0.1)

6 Creating a new variable standardized income

sh$IncZScore<- as.vector(scale(sh$Income, center = TRUE, scale = TRUE))

6.1 creating the new breaks

breaks_Z<- c(-Inf, -2.5,-1, 1, 2.5, Inf)

6.2 creating the custom breaks map

#creating a vector with the colors for the bins in the map
col<- c('#ffffb2','#fecc5c','#fd8d3c','#f03b20','#bd0026')

tm_shape(sh) +
  tm_fill("IncZScore", title= "Standardized Income", breaks=breaks_Z, palette = col)+
  tm_borders()

7 Bivariate Analysis: Analyzing Expenditures by Educational Attainment

A scatter plot and the OLS regression estimates

sh %>% 
  ggplot(aes(x=University, y=Expenses))+
  geom_point()+
  geom_smooth(method = "lm")


lm.sh<-   lm(Expenses ~ University, data=sh)
summary( lm.sh)

Call:
lm(formula = Expenses ~ University, data = sh)

Residuals:
    Min      1Q  Median      3Q     Max 
-132.33  -29.14  -14.76    9.47  293.99 

Coefficients:
            Estimate Std. Error t value          Pr(>|t|)    
(Intercept)    15.69      19.17    0.82              0.42    
University      8.89       0.97    9.16 0.000000000000019 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 64.4 on 88 degrees of freedom
Multiple R-squared:  0.488, Adjusted R-squared:  0.482 
F-statistic: 83.9 on 1 and 88 DF,  p-value: 0.0000000000000193

8 Creating a scatter plot matrix

data<- data.frame( sh$Expenses, sh$University, sh$SecondaryE)
data
chart.Correlation(data, histogram = TRUE, method = "pearson")

8.1 A better option

#pairs.panels(data,
 #            smooth = TRUE,      # If TRUE, draws loess smooths
 #            scale = FALSE,      # If TRUE, scales the correlation text font
 #            density = TRUE,     # If TRUE, adds density plots and histograms
  #           ellipses = TRUE,    # If TRUE, draws ellipses
  #           method = "pearson", # Correlation method (also "spearman" or "kendall")
  #           pch = 21,           # pch symbol
  #           lm = TRUE,         # If TRUE, plots linear fit rather than the LOESS (smoothed) fit
  #           cor = TRUE,         # If TRUE, reports correlations
  #           jiggle = FALSE,     # If TRUE, data points are jittered
  #           factor = 2,         # Jittering factor
  #           hist.col = 4,       # Histograms color
  #           stars = TRUE,       # If TRUE, adds significance level with stars
  #           ci = TRUE)          # If TRUE, adds confidence intervals

That’s all!

LS0tDQp0aXRsZTogIkV4cGxvcmF0b3J5IFNwYXRpYWwgRGF0YSBBbmFseXNpcw0KVG9vbHMgYW5kIFN0YXRpc3RpY3MiDQphdXRob3I6ICJGZWxpcGUgU2FudG9zLU1hcnF1ZXogJiBTZXcgU29vayBZYW4iDQpkYXRlOiAiMTAvMTIvMjAyMCINCm91dHB1dDoNCiAgaHRtbF9ub3RlYm9vazoNCiAgICBjb2RlX2ZvbGRpbmc6IHNob3cNCiAgICBoaWdobGlnaHQ6IG1vbm9jaHJvbWUNCiAgICBudW1iZXJfc2VjdGlvbnM6IHllcw0KICAgIHRoZW1lOiBjb3Ntbw0KICAgIHRvYzogeWVzDQogICAgdG9jX2RlcHRoOiA0DQogICAgdG9jX2Zsb2F0Og0KICAgICAgY29sbGFwc2VkOiBubw0KICAgICAgc21vb3RoX3Njcm9sbDogbm8NCiAgaHRtbF9kb2N1bWVudDoNCiAgICBjb2RlX2Rvd25sb2FkOiB0cnVlDQogICAgZGZfcHJpbnQ6IHBhZ2VkDQogICAgdG9jOiB0cnVlDQogICAgdG9jX2Zsb2F0Og0KICAgICAgY29sbGFwc2VkOiBmYWxzZQ0KICAgICAgc21vb3RoX3Njcm9sbDogZmFsc2UNCiAgICB0b2NfZGVwdGg6IDQNCiAgICBudW1iZXJfc2VjdGlvbnM6IHRydWUNCiAgICBjb2RlX2ZvbGRpbmc6ICJzaG93Ig0KICAgIHRoZW1lOiAiY29zbW8iDQogICAgaGlnaGxpZ2h0OiAibW9ub2Nocm9tZSINCiAgcGRmX2RvY3VtZW50OiBkZWZhdWx0DQogIHdvcmRfZG9jdW1lbnQ6IGRlZmF1bHQNCmJpYmxpb2dyYXBoeTogYmlibGlvLmJpYg0KLS0tDQoNCjxzdHlsZT4NCmgxLnRpdGxlIHtmb250LXNpemU6IDE4cHQ7IGNvbG9yOiBEYXJrQmx1ZTt9IA0KYm9keSwgaDEsIGgyLCBoMywgaDQge2ZvbnQtZmFtaWx5OiAiUGFsYXRpbm8iLCBzZXJpZjt9DQpib2R5IHtmb250LXNpemU6IDEycHQ7fQ0KLyogSGVhZGVycyAqLw0KaDEsaDIsaDMsaDQsaDUsaDZ7Zm9udC1zaXplOiAxNHB0OyBjb2xvcjogIzAwMDA4Qjt9DQpib2R5IHtjb2xvcjogIzMzMzMzMzt9DQphLCBhOmhvdmVyIHtjb2xvcjogIzhCM0E2Mjt9DQpwcmUge2ZvbnQtc2l6ZTogMTJweDt9DQo8L3N0eWxlPg0KDQojIExvYWRpbmcgcmVxdWlyZWQgcGFja2FnZXMNCg0KYGBge3Igc2V0dXB9DQprbml0cjo6b3B0c19jaHVuayRzZXQod2FybmluZyA9IEZBTFNFLCBtZXNzYWdlID0gRkFMU0UsIGVjaG8gPSBGQUxTRSkNCmxpYnJhcnkodGlkeXZlcnNlKQ0KbGlicmFyeShzZikNCmxpYnJhcnkodG1hcCkNCmxpYnJhcnkoUGVyZm9ybWFuY2VBbmFseXRpY3MpICNzY2F0dGVyIHBsb3QgbWF0cml4DQpsaWJyYXJ5KFBlcmZvcm1hbmNlQW5hbHl0aWNzKSAjc2NhdHRlciBwbG90IG1hdHJpeA0KIA0KIyBDaGFuZ2UgdGhlIHByZXNlbnRhdGlvbiBvZiBkZWNpbWFsIG51bWJlcnMgdG8gNCBhbmQgYXZvaWQgc2NpZW50aWZpYyBub3RhdGlvbg0Kb3B0aW9ucyhwcm9tcHQ9IlI+ICIsIGRpZ2l0cz01LCBzY2lwZW49OTk5KQ0KYGBgDQoNCg0KIyBJbXBvcnQgdGhlIHNoYXBlZmlsZQ0KDQoNCmBgYHtyfQ0Kc2g8LSBzdF9yZWFkKCIuLi8uLi9EYXRhL2NpdHkuc2hwIikNCm9iamVjdC5zaXplKHNoKQ0Kcy5zcCA8LSBhcyhzaCwgIlNwYXRpYWwiKQ0KY2xhc3Mocy5zcCkNCmNsYXNzKHNoKQ0KYGBgDQoNCiMjIEluc3BlY3QgdGhlIHRhYmxlIHdpdGhpbiB0aGUgc2hhcGVmaWxlDQoNCmBgYHtyfQ0Kc3RfZHJvcF9nZW9tZXRyeShzaCkNCmBgYA0KDQojIyBBIHF1aWNrIG1hcA0KDQpgYGB7cn0NCnF0bShzaCkNCmBgYA0KIyBDcmVhdGluZyB0aGUgYnJlYWtzIGZvciB0aGUgY3VzdG9tIGJyZWFrcyBtYXAgDQoNCmBgYHtyfQ0KYnJlYWtzPC0gYygtSW5mLCAxNTAwMCwyMDAwMCwyNTAwMCwgSW5mKQ0KYGBgDQoNCiMjIGNyZWF0aW5nIHRoZSBjdXN0b20gYnJlYWtzIG1hcA0KDQpgYGB7cn0NCnRtX3NoYXBlKHNoKSArDQogIHRtX2ZpbGwoIkluY29tZSIsIHRpdGxlPSAiSW5jb21lIixicmVha3M9YnJlYWtzLHBhbGV0dGU9IkdyZWVucyIpKw0KICB0bV9ib3JkZXJzKCkNCmBgYA0KDQoNCiMgQ3JlYXRpbmcgYSBoaXN0b2dyYW0NCg0KYGBge3J9DQpzaCAlPiUgDQogIGdncGxvdChhZXMoeD1JbmNvbWUpKSsNCiAgZ2VvbV9oaXN0b2dyYW0oYmlucz0xMCwgY29sb3I9ImJsYWNrIiwgZmlsbD0id2hpdGUiKQ0KYGBgDQoNCg0KIyBjcmVhdGluZyBhIGJveCBwbG90DQoNCmdlb21faml0dGVyKHNpemUgPSAzLCBhbHBoYSA9IDAuMywgd2lkdGggPSAwLjEpIGlzIHVzZWQgc28gdGhhdCB0aGUgcG9pbnRzIGFyZSBub3QgYWxsIGNyYW1tZWQgaW4gdGhlDQpzYW1lIGxpbmUNCg0KYGBge3J9DQpzZXQuc2VlZCgyMDE5KQ0Kc2ggJT4lIA0KICBnZ3Bsb3QoYWVzKHg9SW5jb21lLCB5PU11bmljaXBhbGkpKSsNCiAgZ2VvbV9ib3hwbG90KCkrDQogIGdlb21faml0dGVyKHNpemUgPSAzLCBhbHBoYSA9IDAuMywgd2lkdGggPSAwLjEpDQpgYGANCg0KIyBDcmVhdGluZyBhIG5ldyB2YXJpYWJsZSBzdGFuZGFyZGl6ZWQgaW5jb21lDQoNCmBgYHtyfQ0Kc2gkSW5jWlNjb3JlPC0gYXMudmVjdG9yKHNjYWxlKHNoJEluY29tZSwgY2VudGVyID0gVFJVRSwgc2NhbGUgPSBUUlVFKSkNCmBgYA0KDQojIyBjcmVhdGluZyB0aGUgbmV3IGJyZWFrcw0KDQpgYGB7cn0NCmJyZWFrc19aPC0gYygtSW5mLCAtMi41LC0xLCAxLCAyLjUsIEluZikNCmBgYA0KDQojIyBjcmVhdGluZyB0aGUgY3VzdG9tIGJyZWFrcyBtYXANCg0KYGBge3J9DQojY3JlYXRpbmcgYSB2ZWN0b3Igd2l0aCB0aGUgY29sb3JzIGZvciB0aGUgYmlucyBpbiB0aGUgbWFwDQpjb2w8LSBjKCcjZmZmZmIyJywnI2ZlY2M1YycsJyNmZDhkM2MnLCcjZjAzYjIwJywnI2JkMDAyNicpDQoNCnRtX3NoYXBlKHNoKSArDQogIHRtX2ZpbGwoIkluY1pTY29yZSIsIHRpdGxlPSAiU3RhbmRhcmRpemVkIEluY29tZSIsIGJyZWFrcz1icmVha3NfWiwgcGFsZXR0ZSA9IGNvbCkrDQogIHRtX2JvcmRlcnMoKQ0KYGBgDQoNCiMgQml2YXJpYXRlIEFuYWx5c2lzOiBBbmFseXppbmcgRXhwZW5kaXR1cmVzIGJ5IEVkdWNhdGlvbmFsIEF0dGFpbm1lbnQgDQoNCkEgc2NhdHRlciBwbG90IGFuZCB0aGUgT0xTIHJlZ3Jlc3Npb24gZXN0aW1hdGVzICANCg0KYGBge3J9DQpzaCAlPiUgDQogIGdncGxvdChhZXMoeD1Vbml2ZXJzaXR5LCB5PUV4cGVuc2VzKSkrDQogIGdlb21fcG9pbnQoKSsNCiAgZ2VvbV9zbW9vdGgobWV0aG9kID0gImxtIikNCg0KbG0uc2g8LSAgIGxtKEV4cGVuc2VzIH4gVW5pdmVyc2l0eSwgZGF0YT1zaCkNCnN1bW1hcnkoIGxtLnNoKQ0KYGBgDQoNCiMgQ3JlYXRpbmcgYSBzY2F0dGVyIHBsb3QgbWF0cml4DQoNCmBgYHtyfQ0KZGF0YTwtIGRhdGEuZnJhbWUoIHNoJEV4cGVuc2VzLCBzaCRVbml2ZXJzaXR5LCBzaCRTZWNvbmRhcnlFKQ0KZGF0YQ0KY2hhcnQuQ29ycmVsYXRpb24oZGF0YSwgaGlzdG9ncmFtID0gVFJVRSwgbWV0aG9kID0gInBlYXJzb24iKQ0KYGBgDQoNCiMjIEEgYmV0dGVyIG9wdGlvbiANCg0KYGBge3J9DQojcGFpcnMucGFuZWxzKGRhdGEsDQogIyAgICAgICAgICAgIHNtb290aCA9IFRSVUUsICAgICAgIyBJZiBUUlVFLCBkcmF3cyBsb2VzcyBzbW9vdGhzDQogIyAgICAgICAgICAgIHNjYWxlID0gRkFMU0UsICAgICAgIyBJZiBUUlVFLCBzY2FsZXMgdGhlIGNvcnJlbGF0aW9uIHRleHQgZm9udA0KICMgICAgICAgICAgICBkZW5zaXR5ID0gVFJVRSwgICAgICMgSWYgVFJVRSwgYWRkcyBkZW5zaXR5IHBsb3RzIGFuZCBoaXN0b2dyYW1zDQogICMgICAgICAgICAgIGVsbGlwc2VzID0gVFJVRSwgICAgIyBJZiBUUlVFLCBkcmF3cyBlbGxpcHNlcw0KICAjICAgICAgICAgICBtZXRob2QgPSAicGVhcnNvbiIsICMgQ29ycmVsYXRpb24gbWV0aG9kIChhbHNvICJzcGVhcm1hbiIgb3IgImtlbmRhbGwiKQ0KICAjICAgICAgICAgICBwY2ggPSAyMSwgICAgICAgICAgICMgcGNoIHN5bWJvbA0KICAjICAgICAgICAgICBsbSA9IFRSVUUsICAgICAgICAgIyBJZiBUUlVFLCBwbG90cyBsaW5lYXIgZml0IHJhdGhlciB0aGFuIHRoZSBMT0VTUyAoc21vb3RoZWQpIGZpdA0KICAjICAgICAgICAgICBjb3IgPSBUUlVFLCAgICAgICAgICMgSWYgVFJVRSwgcmVwb3J0cyBjb3JyZWxhdGlvbnMNCiAgIyAgICAgICAgICAgamlnZ2xlID0gRkFMU0UsICAgICAjIElmIFRSVUUsIGRhdGEgcG9pbnRzIGFyZSBqaXR0ZXJlZA0KICAjICAgICAgICAgICBmYWN0b3IgPSAyLCAgICAgICAgICMgSml0dGVyaW5nIGZhY3Rvcg0KICAjICAgICAgICAgICBoaXN0LmNvbCA9IDQsICAgICAgICMgSGlzdG9ncmFtcyBjb2xvcg0KICAjICAgICAgICAgICBzdGFycyA9IFRSVUUsICAgICAgICMgSWYgVFJVRSwgYWRkcyBzaWduaWZpY2FuY2UgbGV2ZWwgd2l0aCBzdGFycw0KICAjICAgICAgICAgICBjaSA9IFRSVUUpICAgICAgICAgICMgSWYgVFJVRSwgYWRkcyBjb25maWRlbmNlIGludGVydmFscw0KDQoNCmBgYA0KDQpUaGF0J3MgYWxsIQ0KDQohW10oaHR0cHM6Ly9tZWRpYS5naXBoeS5jb20vbWVkaWEveFVPeGY3WGZtcHh1U29kZTFPL2dpcGh5LmdpZikNCg0KDQoNCg0KDQo=